function best_mix = clusteringCEM(mix, dataset, options)

% initialisation
[N D] = size(dataset.Y);
switch mix.covar_type
    case 'full'
        npars = D + (D*(D+1)/2);
    case 'diagonal'
        npars = D + D;
    case 'spherical'
        npars = D + 1;
    otherwise
        error(['Unknown covariance type ', covar_type]);
end
nparsover2 = npars/2;
lmin = realmax;
best_mix = mix;


% first expectation step
[post act] = gmmpostRobust(mix, dataset.Y);

if options.verbose
    fprintf('initial negloglikelihood = %10.10f\n', -estimateLoglikelihood(mix, act));
    fprintf('initial negloglikelihood + descriptionlength = %10.10f\n', -estimateLoglikelihood(mix, act)+computeDescriptionLength(mix, nparsover2, N));
end


% the outer loop
k_cont = 1;
while k_cont
    % init negloglikelihoods
    iters_count = 0;
    curr_negloglikelihood = realmax/2;
    prev_negloglikelihood = realmax;

    % the inner loop
    while abs(prev_negloglikelihood-curr_negloglikelihood) > abs(options.th*prev_negloglikelihood)
        iters_count = iters_count + 1;
        prev_negloglikelihood = curr_negloglikelihood;
        comps_to_kill = zeros(1, mix.ncentres);

        % the inner inner loop
        for k = 1:mix.ncentres
            % compute posterior probability
            valid_inds = find(comps_to_kill==0);
            post(:, valid_inds) = act(:, valid_inds).*repmat(mix.priors(valid_inds), N, 1);

            % normalise posterior probability
            post_sum = sum(post(:, valid_inds), 2);
            post_sum = post_sum + (post_sum==0);
            post(:, post_sum==0) = 1/mix.ncentres;
            post(:, k) = post(:, k)./post_sum;

            % update mixing probabilities
            if isfield(options, 'no_kill') && options.no_kill == 1
                mix.priors(k) = max(sum(post(:, k)), 0)/N;
            else
                mix.priors(k) = max(sum(post(:, k))-nparsover2, 0)/N;
            end
            mix.priors = mix.priors/sum(mix.priors);

            % either set component to be killed or do maxmisation
            if mix.priors(k) == 0
                if options.verbose
                    fprintf('X');
                end

                comps_to_kill(k) = 1;
            else
                if options.verbose
                    fprintf('*');
                end

                % do maxmisation
                [mu sigma] = computeMLEGaussian(dataset, post(:, k), mix.covar_type, options.reg*0.1);
                mix.centres(k, :) = mu;
                mix.covars(:, :, k) = sigma;

                % update class conditional probability
                act(:, k) = gaussRobust(mix.centres(k, :), mix.covars(:, :, k), dataset.Y);
            end
        end

        % kill components
        temp_vars.a = 0;
        mix = killComponents(mix, temp_vars, find(comps_to_kill==1));
        act = act(:, comps_to_kill==0);
        post = act.*repmat(mix.priors, N, 1);

        curr_negloglikelihood = -estimateLoglikelihood(mix, act);
        dl = computeDescriptionLength(mix, nparsover2, N);
        if options.verbose
            tol = (prev_negloglikelihood-curr_negloglikelihood) / abs(options.th*prev_negloglikelihood);
            fprintf('ITERS (%d)   (%f)   negloglikelihood: %10.10f   (+ %10.10f = %10.10f)   knz: %d \n', ...
                     iters_count, tol, curr_negloglikelihood, dl, curr_negloglikelihood+dl, mix.ncentres);
        end
    end

    if isfield(options, 'no_mml') && options.no_mml == 1
        curr_loglikelihood = -estimateLoglikelihood(mix, act);
    else
        curr_loglikelihood = -estimateLoglikelihood(mix, act) + computeDescriptionLength(mix, nparsover2, N);
    end
    if curr_loglikelihood < lmin
        lmin = curr_loglikelihood;
        best_mix = mix;
    end

    % if there are more than 1 component, kill least significant component
    if mix.ncentres > options.kmin
        [temp least_comp] = min(mix.priors);
        keep_comps = setdiff(1:mix.ncentres, least_comp);

        temp_vars.a = 0;
        mix = killComponents(mix, temp_vars, least_comp);
        post = post(:, keep_comps);
        act = act(:, keep_comps);
    else
        k_cont = 0;
    end

    if options.verbose
        fprintf('%10.10f   %d\n', lmin, best_mix.ncentres);
    end
end
